function out = getImpulseResponses(outEstim,IRFlength,shockSizeScalar,xfGIRF,xsGIRF,plotGIRFsON)
if nargin == 5
    plotGIRFsON = 1;
end


%% Impulse response functions
model       = outEstim.model;
ne          = size(model.eta,2);
g0          = model.g0;
gx          = model.gx;
gxx         = model.gxx;
labely      = model.labely;
% Add short-yield
nn = model.ny + 1;
labely    = [labely ,{'$r_t$'}];
g0(nn,1)  = model.by0(1,1);
gx(nn,:)  = model.byx(1,:);
gxx(nn,:,:) = model.byxx(1,:,:);
% Add long yield
maxMat = max(outEstim.model.maturities);
nn = nn + 1;
labely     = [labely ,{['$r^{(',num2str(maxMat),')}_t$']}];
g0(nn,1)  = model.by0(maxMat,1);
gx(nn,:)  = model.byx(maxMat,:);
gxx(nn,:,:) = model.byxx(maxMat,:,:);
% Add slope 
nn = nn + 1;
labely     = [labely ,{['$r^{(',num2str(maxMat),')}_t-r^{(4)}_t$']}];
g0(nn,1)  = model.by0(maxMat,1)-model.by0(4,1);
gx(nn,:)  = model.byx(maxMat,:)-model.byx(4,:);
gxx(nn,:,:) = model.byxx(maxMat,:,:)-model.byxx(4,:,:);

% Conditional expectations of bond yields at max maturity for m periods ahead
m      = outEstim.setupFilter.CSregress_m;
maxMat = max(outEstim.model.maturities);
var0   = model.by0(maxMat-m,1);
varx   = model.byx(maxMat-m,:);
varxx  = model.byxx(maxMat-m,:,:);
varxxx = model.byxxx(maxMat-m,:,:,:);
[Er0,Erx,Erxx]= projectionStepExpectedControlClosedForm(var0,varx,varxx,varxxx,...
       m,model,outEstim.setupFilter.setupProj);
nn = nn + 1;
labely    = [labely ,{['$E_{t}\left[r^{(',num2str(maxMat-m),')}_{t+',num2str(m),'}\right]','-r^{(',num2str(maxMat),')}_t$']}];
g0(nn,1)  = Er0(end,1)-model.by0(maxMat,1);
gx(nn,:)  = Erx(end,:)-model.byx(maxMat,:);
gxx(nn,:,:) = Erxx(end,:,:)-model.byxx(maxMat,:,:);

% Adding term premium
nn = nn + 1;
labely    = [labely ,{['$TP^{(',num2str(maxMat),')}_t$']}];
g0(nn,1)  = model.TP0(maxMat,1);
gx(nn,:)  = model.TPx(maxMat,:);
gxx(nn,:,:) = reshape(model.TPxx(maxMat,:,:),1,model.nx,model.nx);

%adding sdf
labely    = [labely ,{'$M_{t,t+1}$'},{'$M^{noEZ}_{t,t+1}$'},{'$M^{noD}_{t,t+1}$'}];

% Test the size of xfGIRF, the point around which the GIRFs are
% computed
nx  = size(model.hx,1);
numXf = 1;
if ~isempty(xfGIRF)
    if size(xfGIRF,1) ~= nx
        error('xfGIRF must be of size nx*1')
    end
    if size(xfGIRF,2) > 1
        % we have many points for xf to evaluate the GIRFs at
        numXf = size(xfGIRF,2);
        IRFySave = zeros(size(gx,1)+3,IRFlength,ne,numXf);
    else
        % We have only one point to evaluate the GIRFs at
        IRFySave = zeros(size(gx,1)+3,IRFlength,ne,1);
    end
else
    IRFySave = zeros(size(gx,1)+3,IRFlength,ne,1);
end
parfor j=1:ne
    j
    shockSize      = zeros(ne,1);
    shockSize(j,1) = shockSizeScalar;
    if numXf == 1
        for s=1
            % Impulse responses for the SDF
            [IRFsdf,IRFsdfnoEZ,IRFsdfnoD] = computeIRFsdf(model,shockSize,xfGIRF,xsGIRF,model.appOrder,IRFlength);

            [IRFy,IRFx] = GIRFPruning_Proj(gx,gxx,model.h0,model.hx,model.hxx,model.eta,...
                shockSize,xfGIRF,model.appOrder,IRFlength);
            IRFySave(:,:,j,s) = [IRFy(:,:,model.appOrder); IRFsdf';IRFsdfnoEZ';IRFsdfnoD'];
        end
    else
        for s=1:numXf
            % Impulse responses for the SDF
            [IRFsdf,IRFsdfnoEZ,IRFsdfnoD] = computeIRFsdf(model,shockSize,xfGIRF(:,s),xsGIRF(:,s),model.appOrder,IRFlength);

            [IRFy,IRFx] = GIRFPruning_Proj(gx,gxx,model.h0,model.hx,model.hxx,model.eta,...
                shockSize,xfGIRF(:,s),model.appOrder,IRFlength);
            IRFySave(:,:,j,s) = [IRFy(:,:,model.appOrder); IRFsdf';IRFsdfnoEZ';IRFsdfnoD'];
        end
    end

    % Extra scaling of y
    IRFy            = [IRFy;repmat(IRFsdf',1,1,2);repmat(IRFsdfnoEZ',1,1,2);repmat(IRFsdfnoD',1,1,2)];
    scalingY = ones(size(IRFy,1),1);
    % Showing only first element in x
    if plotGIRFsON == 1 && numXf == 1
        IRFx(2:end,:,:) = IRFx(2:end,:,:)*0;
        plotImpulseResponse(repmat(scalingY,1,size(IRFy,2),size(IRFy,3)).*IRFy*100,...
            IRFx*100,labely,model.labelx,shockSize,model.appOrder);
    end

end

% saving output
if numXf == 1
    out.IRFy  = squeeze(IRFySave);
else
    out.IRFy  = IRFySave;
end
out.labely = labely;
out.xfGIRF = xfGIRF;
out.shockSizeScalar = shockSizeScalar;
